\documentclass[a4paper,10pt]{article}
\usepackage[pdftex]{graphicx}
\usepackage{polski}
\usepackage[utf8]{inputenc}
\usepackage{textcomp}
\usepackage{pifont}
\usepackage{color}
\usepackage{listings}
\usepackage{amsfonts}

% opcje listingów
\definecolor{darkgray}{rgb}{0.95,0.95,0.95}
\lstset{language=c++}
\lstset{backgroundcolor=\color{darkgray}}

\title{\textbf{Elementy Bioinformatyki \\ Projekt 2\\ \vspace{2 mm} {\large Zadanie 7. Sekwencjonowanie przez hybrydyzację}}}
\author{\textbf{Przemysław Chmielewski} [125543]\\ \textbf{Radosław Czubala} [125552]\\ \textbf{Przemysław Konopa} [1199]}
\date{23 stycznia 2014r.}

\begin{document}
 
\maketitle

\section{Treść zadania}
,,Dany jest multizbiór słów tej samej 
długości k nad alfabetem DNA, sprawdź czy istnieje (i w tym przypadku zwróć na wyjście) 
sekwencja, dla których podane wejście stanowiło spektrum. Co więcej, program powinien być 
odporny na podaną przez użytkownika w zakresie 0-3 liczbę błędów (pozytywnych i 
negatywnych) i mieć (przy ustalonym k) złożoność wielomianową.''
\section{Wprowadzenie do zagadnienia}

Problem sekwencjonowania łańcuchów DNA jest w ogólności problemem silnie NP-trudnym. W problemie tym mamy
dany zbiór (spektrum) $\mathbb{S}$ słów (oligonukleotydów) o jednakowej długości $l$ nad alfabetem $\{A, C, G, T\}$.
Ponadto mamy podaną długość sekwencji oryginalnej $n$. Celem sekwencjonowania jest odtworzenie oryginalnego łańcucha DNA, 
który został poddany hybrydyzacji, na podstawie powyżej zdefiniowanych danych wejściowych.
Dla ułatwienia dalszego opisu wprowadźmy kilka terminów: 
\begin{description}
  \item[Oligonukleotyd] fragment kwasu nukleinowego długości od kilku do kilkunastu nukleotydów (słowo o długości $l$ nad alfabetem $\{A, C, G, T\}$)
  \item[Spektrum] zbiór $\mathbb{S}$ oligonukleotydów
  \item[Sekwencja DNA] słowo o długości $n$ nad alfabetem $\{A, C, G, T\}$
  \item[Błąd negatywny] występuje gdy w spektrum nie ma niektórych oligonukleotydów
  \item[Błąd pozytywny] występuje gdy w spektrum znajdują się dodatkowe oligonukleotydy
  \item[Odległość między oligonukleotydami] definiuje się między każdą parą słów jako najmniejsze przesunięcie jednego oligonukletydu
                                            względem drugiego, które umożliwi połączenie ich w jeden łańcuch. Dla przykładu słowo $ACGTA$
                                            znajduje się w odległości 2 względem słowa $ACACG$ (relacja w drugą stronę - słowo $ACACG$ 
                                            znajduje się w odległości 4 względem słowa $ACGTA$).

  \begin{figure}[h]
    \footnotesize\centering
    \includegraphics[width=0.8\textwidth,keepaspectratio]{Distance1.png}
  \end{figure}

  \begin{figure}[h]
    \footnotesize\centering
    \includegraphics[width=0.8\textwidth,keepaspectratio]{Distance2.png}
  \end{figure}
  \item[Graf DNA] skierowany graf w którym wierzchołkami są oligonukleotydy, a krawędzie odpowiadają odległością między nimi (zdefiniowana powyżej)
\end{description}

W zaproponowanych przez nas podejściach będziemy wykorzystywać teorię grafów w celu szybkiego łączenia oligonukleotydów oddalonych od siebie o jak najmniejszą wartość.

W analizie problemu można wyróżnić cztery szczególne przypadki - idealny, z błędami negatywnymi, z błędami pozytywnymi, z obydwoma typami błędów.
Przypadek idealny to taki, gdzie na mikromacierzy dopasowały się wszystkie słowa występujące w sekwencji.
Mówimy, że dane odwzorowanie łańcucha DNA na mikromacierzy zawiera błędy negatywne, jeśli z jakiegoś powodu niektóre słowa nie
zostały dopasowane do macierzy. Otrzymujemy więc niepełne spektrum sekwencji.
Błędami pozytywnymi określamy nadmiar informacji, tzn. na mikromacierzy zostały dopasowane słowa, które nie występują w sekwencji
oryginalnej. Mamy więc zbiór słów które są nadzbiorem zbioru słów wchodzących w skład oryginalnej sekwencji.

Przy projektowaniu heurystyk do problemu SBH założyliśmy że $n << 4^l$.

\section{Generowanie losowych instancji}
Wszystkie algorytmy heurystyczne były sprawdzane na instancjach testowych. 
Aby jednak nie ograniczać się do zamkniętej puli testów, postanowiliśmy 
napisać generator losowych instancji. 

Generator instacji pozwala na wygenerowanie instacji tylko z błędami pozytywnymi
lub tylko z błędami negatywnymi. W obu przypadkach parametry to:
\begin{itemize}
 \item wielkość spektrum idealnego $|S|$
 \item ilość błędów $e$
 \item długość oligonukleotydu $l$
\end{itemize}

Pierwszym krokiem wykorzystywanym przy generowaniu obu typów instacji jest stworzenie
losowego łańcucha DNA o długości $k$, który nie zawiera błędów negatywnych (żaden 
podciąg o długości $l$ się nie powtarza). Najprostszym sposobem jest wygenerowanie 
ciągu o długości $k$ składającego się z losowych liter A, C, T, G. Ponieważ możemy 
założyć że $k$ będzie dużo mniejsze niż $4^l$ prawdopodobieństwo że tak wygenerowany
ciąg będzie niepoprawny jest bardzo małe. Gdyby jednak okazało się że błąd wystąpił,
należy uruchomić metodę od początku.

Aby wygenerować ciąg którego spektrum idealne ma rozmiar $|S|$ z $e$ błędami pozytywnymi
wystarczy wygenerować łańcuch DNA o długości $|S|+l-1$ i do tak utworzonego spektrum 
dodać $e$ unikalnych oligonukleotydów.

Wygenerowanie instancji z błędami negatywnymi jest jeszcze prostsze. Najpierw należy
wygenerować łańcuch DNA o długości $|S|+e+l-1$, a następnie usunąć z tego spektrum
$e$ losowych oligonukleotydów.

\section{SBH z błędami pozytywnymi}
Jak już wcześniej mówiliśmy, sekwencjonowanie łańcucha DNA z błędami pozytywnymi charakteryzuje się
nadmiarową ilością słów otrzymanych w spektrum. Mamy więc pewność, że pośród tych słów znajdzie się dokładnie $n-l+1$ z których będziemy mogli zrekonstruować oryginalną sekwencję.

Pierwszym krokiem jest zbudowanie grafu DNA $G$ który zawiera tylko łuki o wadze 1. Graf taki posiada dokładnie $|S|$ wierzchołków a ilość krawędzi jest mniejsza niż $4 \cdot |S|$.
Gdy graf $G$ jest grafem acyklicznym, wtedy używając sortowania topologicznego oraz programowania dynamicznego możemy dla każdego wierzchołka $u$ wyliczyć 
najdłuższą ścieżkę zaczynającą się w tym wierzchołku - oznaczmy ją przez $distance[u]$. Taką tablicę można wyliczyć w $O(|S|)$. Używając tej tablicy 
możemy znaleźć najdłuższą ścieżkę w grafie $G$ która jest zarazem optymalnym rozwiązaniem dla problemu SBH z błędami pozytywnymi.

Gdy graf nie jest acykliczny, nasza heurystyka usuwa krawędzie używając poniższego algorytmu:

Algorytm 1
\begin{enumerate}
 \item Znajdź silnie spójne składowe grafu $G$
 \item Jeśli graf jest acykliczny, zakończ algorytm
 \item Stwórz nowy graf $G_{SCC}$, traktując każdą silnie spójną składową $G$ jako wierzchołek oraz
       tworząc łuk z silnie spójnej składowej $SCC_1$ do $SCC_2$ jeśli jakikolwiek wierzchołek z $SCC_1$ łączy się z dowolnym wierzchołkiem z $SCC_2$
 \item Posortuj graf $G_{SCC}$ topologicznie
 \item Znajdź pierwszą silnie spójną składową $SCC$ (w kolejności sortowania topologicznego) która zawiera więcej niż jeden wierzchołek
 \item Dla danej silnie spójnej składowej $SCC$ znajdź należący do niej wierzchołek $k$, który łączy się z wierzchołkiem $v$ z 
       innej silnie spójnej składowej o największej wartości $distance[v]$
 \item Usuń wszystkie połączenia wychodzące z $k$, które łączą się z innymi wierzchołkami wchodzącymi w skład $SCC$
 \item Wróć do 1
\end{enumerate}

Wykonanie jednej iteracji tego algorytmu (aż do punktu 8) wymaga $O(|S|)$ operacji. W każdej
iteracji usuwana jest przynajmniej jedna krawędź w grafie, a więc maksymalnie wykonamy $O(|S|)$ iteracji. Stąd ostateczna złożoność powyższego
algorytmu to $O(|S|^2)$.

Mając do dyspozycji algorytm usuwania cykli w grafie proponujemy następujący algorytm rozwiązujący problem sekwencjonowania DNA z błędami pozytywnymi:

Algorytm 2
\begin{enumerate}
 \item Utwórz graf skierowany $G$. Niech zbiór wierzchołków $V$ stanowią słowa tworzące spektrum sekwencji DNA. Dla każdej pary 
 wierzchołków $(i,j)$ dodaj łuk od wierzchołka $i$ do $j$ $\iff$ gdy słowo w wierzchołku $j$ jest oddalone o 1 od słowa w wierzchołku $i$
 \item Usuń cykle w grafie $G$ tworząc graf $G_{DAG}$ korzystając z powyższego algorytmu
 \item Uszereguj topologicznie $G_{DAG}$
 \item Wybierz najdłuższą możliwą ścieżkę w grafie $G_{DAG}$
\end{enumerate}

Złożoność poszczególnych kroków algorytmu jest następująca:
\begin{itemize}
 \item $O(|S|^2)$ dla punktu 1
 \item $O(|S|^2)$ dla punktu 2
 \item $O(|S|)$ dla punktu 3
 \item $O(|S|)$ dla punktu 4
\end{itemize}

Ostateczna złożoność całej heurystyki wynosi $O(|S|^2)$.

\section{SBH dla przypadku ogólnego}
W przypadku podejścia ogólnego nie mamy gwarancji, że istnieje ścieżka po łukach o koszcie 1, która utworzy sekwencję odpowiedniej długości. 
Musimy więc uwzględnić mniej korzystne połączenia o kosztach 2 i wyższych. W związku z tym wykorzystamy algorytm sekwencjonowania z błędami 
pozytywnymi uogólniając jego działanie. Algorytm po modyfikacji wygląda następująco:

Algorytm 3
\begin{enumerate}
 \item Wczytaj zbiór wszystkich słów tworzących spektrum do słownika \texttt{OLIGS}
 \item Utwórz graf skierowany $G$. Niech zbiór wierzchołków $V$ stanowią słowa w słowniku \texttt{OLIGS}. Dla każdej pary wierzchołków $(i,j)$ dodaj 
       łuk od wierzchołka $i$ do $j$ o koszcie równym odległości słowa w wierzchołku $j$ względem słowa w wierzchołku $i$
 \item Znajdź najniższy koszt łuku \texttt{MIN} w grafie $G$
 \item Utwórz graf skierowany $G'$, który jest podgrafem grafu $G$ i zawiera tylko łuki o koszcie \texttt{MIN}
 \item Usuń cykle w grafie $G'$ tworząc graf $G'_{DAG}$ korzystając z Algorytmu 1
 \item Uszereguj topologicznie $G'_{DAG}$
 \item Wybierz najdłuższą możliwą ścieżkę \texttt{PATH} w grafie $G'_{DAG}$
 \item Usuń ze słownika \texttt{OLIGS} słowa wchodzące w skład ścieżki \texttt{PATH}, w zamian dodając słowo utworzone przez połączenie słów zgodnie z kolejnością w ścieżce \texttt{PATH}
 \item Jeśli słownik \texttt{OLIGS} zawiera dokładnie jedno słowo, to zakończ, w przeciwnym wypadku idź do punktu 2
\end{enumerate}

Z punktu widzenia złożoności należy zauważyć, że punkt 3 można zrealizować w ramach punktu 2. Zauważmy ponadto, że punkty od 4 do 7 to dokładnie algorytm dla przypadku błędów pozytywnych, 
z uwzględnieniem łuków o koszcie MIN (zamiast o koszcie 1 jak podano w poprzednim algorytmie). Złożoności pozostałych punktów to:
\begin{itemize}
 \item $O(S^2)$ dla punktu 2
 \item $O(1)$ dla usuwania w punkcie 8, tworzenie słowa jest liniowo zależne od ilości oligonukleotydów, które wchodzą w jego skład
\end{itemize}
W związku z tym ostateczna złożoność algorytmu to $O(|S|^3)$. Przy czym $|S|^3$ występuje w najgorszym wypadku, 
kiedy w grafie wyjściowym będziemy mieć słowa o długości równej liczbie słów instancji i będą istnieć tylko pojedyncze połączenia o kosztach $1, 2, 3, ..., S$.

Ideą algorytmu jest wyszukanie teoretycznie najkorzystniejszych połączeń o koszcie 1 tworzących sekwencje 
pełne (tj. w ich skład wchodzą tylko istniejące słowa), a jeśli nie uda się utworzyć na ich podstawie sekwencji wyjściowej, 
to połączenie sekwencji utworzonych w ten sposób uwzględniając brakujące słowa, czyli łącząc łukami o koszcie 2 i wyższymi. 

Przykład:
$$OLIGS = { AGC, GCT, CTG, TGT, GTG, TGA, GAC, ACT }$$
$$n = 10$$

\begin{figure}[h]
  \footnotesize\centering
  \includegraphics[width=\textwidth,keepaspectratio]{Graph1.png}
\end{figure}

Powyżej zaprezentowano graf utworzony w kroku 4 dla wartości $MIN=1$. Wcześniejsze kroki są jednoznacznie interpretowalne, 
a utworzenie pełnego grafu z punktu 2 bardzo zamgliłoby obraz grafu. Usuńmy teraz cykle.
Jak widać w powyższym grafie mamy następujące silnie spójne składowe:
\begin{itemize}
 \item $A$ - {AGC}
 \item $B$ - {GCT}
 \item $C$ - {CTG, TGT, GTG, TGA, GAC, ACT}
\end{itemize}
Musimy więc rozbić wszystkie silnie spójne składowe usuwając cykle, co umożliwi nam utworzenie najdłuższej możliwej ścieżki od ustalonego punktu.
Uzyskaliśmy tylko jedną silnie spójną składową $C$, która posiada więcej niż jeden wierzchołek, więc ją należy rozbić. 
Zauważmy, że po uszeregowaniu topologicznie spójnych składowych uzyskamy następujące uporządkowanie: $A \rightarrow B \rightarrow C$, 
a więc bierzemy pierwszy lepszy wierzchołek silnie spójnej składowej C i usuwamy z niego wszystkie łuki. 
Ponieważ wybieramy dowolny wierzchołek, załóżmy, że jest to wierzchołek $TGT$, usuwamy zatem łuk prowadzący do $GTG$.
Ponownie generujemy silnie spójne składowe grafu:
\begin{itemize}
 \item $A$ - {AGC}
 \item $B$ - {GCT}
 \item $C$ - {CTG, TGA, GAC, ACT}
 \item $D$ - {TGT}
 \item $E$ - {GTG}
\end{itemize}
Ponieważ znów mamy silnie spójną składową, która ma więcej niż jeden wierzchołek, to z niej będziemy usuwać cykl. 
Tym razem jednak po uszeregowaniu topologicznie silnie spójnych składowych otrzymujemy zależność $A \rightarrow B \rightarrow C \rightarrow D \rightarrow E$. 
Oznacza to, że w silnie spójnej składowej C musimy odnaleźć wierzchołek prowadzący do silnie spójnej składowej $D$. 
Takim wierzchołkiem jest $CTG$, usuwamy więc z niego wszystkie łuki wychodzące z niego, które prowadzą do innych 
wierzchołków silnie spójnej składowej $C$. Zostanie usunięty łuk $CTG \rightarrow TGA$.
W kolejnej iteracji podziału na silnie spójne składowe uzyskamy siedem jedno-elementowych zbiorów, a więc usuwanie 
cykli możemy zakończyć. Otrzymujemy graf pokazany na następnej stronie.

\begin{figure}[h]
  \footnotesize\centering
  \includegraphics[width=\textwidth,keepaspectratio]{Graph2.png}
\end{figure}

Po uszeregowaniu topologicznie możemy uzyskać m.in. następującą kolejność wierzchołków:
$AGC, GCT, GTG, TGA, GAC, ACT, CTG, TGT$

Szukamy najdłuższej możliwej ścieżki korzystając z uporządkowania topologicznego: $GTG, TGA, GAC, ACT, CTG, TGT$
i tworzymy nowe słowo sklejając poszczególne elementy ze sobą w sekwencję $GTGACTGT$.
Usuwamy wierzchołki, z których utworzyliśmy powyższą sekwencją. Dodajemy otrzymaną sekwencję do grafu, który obecnie wygląda następująco:

\begin{figure}[h]
  \footnotesize\centering
  \includegraphics[width=\textwidth,keepaspectratio]{Graph3.png}
\end{figure}

Uwzględnione są tylko łuki o koszcie $MIN=1$ (czyli niejako uzyskaliśmy las grafów). 
Ponieważ graf jest acykliczny to możemy przystąpić do łączenia elementów. 
Najdłuższą ścieżką jest $AGC, GCT$ (liczy się ilość wierzchołków wchodzących w skład ścieżki, nie łączna długość słów!), 
zatem tą sekwencję łączymy i dodajemy do grafu usuwając wykorzystane do jej utworzenia wierzchołki.
Otrzymujemy wówczas poniższy graf:

\begin{figure}[h]
  \footnotesize\centering
  \includegraphics[width=\textwidth,keepaspectratio]{Graph4.png}
\end{figure}
		 
Jest to graf pełny, nad łukami zaznaczono ich koszty. Ponieważ najniższym kosztem jest $MIN=4$, oznacza to, że nie istnieje już 
żadne bezpośrednie połączenie, musimy połączyć te sekwencje "grubymi nićmi" aby utworzyć sekwencję końcową:
$AGCTGTGACTGT$
Zauważmy, że jest ona dłuższa, niż $n=10$. W związku z tym sekwencję wynikową przycinamy do długości odpowiadającej spodziewanemu wynikowi otrzymując w rezultacie:
$AGCTGTGACT$.

Algorytm zakończył swoje działanie. Zupełnie przypadkowo uzyskana sekwencja zgadza się z sekwencją wyjściową, mimo tego, że w ostatnim kroku algorytmu zostały sklejone słowa w sposób niekorzystny.
Należy zauważyć, że znaczny wpływ na działanie algorytmu ma sposób wyboru wierzchołka w silnie spójnej składowej grafu z
którego zostaną usunięte łuki, a także kolejność występowania wierzchołków w silnie spójnych składowych grafu. Jeśli na 
początku zamiast usuwać łuk z wierzchołka $TGT$ usunęlibyśmy dowolny inny łuk z tej silnie spójnej składowej grafu, to dalsze iteracje miałyby inny przebieg.

\section{SBH dla błędów negatywnych}
W algorytmie dla błędów pozytywnych korzystaliśmy z funkcji odległości dwóch łańcuchów tekstowych względem siebie. Było 
to możliwe dzięki temu, że porównywane łańcuchy miały tą samą długość. W przypadku algorytmu ogólnego, w którym również 
wykorzystujemy pojęcie odległości do konstrukcji łuków wkrada się pewna nieścisłość. Otóż porównywane są łańcuchy o 
różnych długościach. Jeśli więc mamy dwa łańcuchy $S_1$=$AAAAAAA$ (długości 7 znaków) oraz $S_2$=$GGGGGGGGAAA$ (długości 11 znaków), 
to zostanie utworzony łuk od $S_1$ do $S_2$ o koszcie 7, a w drugą stronę od $S_2$ do $S_1$ o koszcie 8. Zgodnie z algorytmem 
ogólnym uwzględniony zostanie łuk o niższym koszcie jako lepszy. Jak widać jest to złudne wrażenie, ponieważ łącząc elementy 
przy pomocy drugiego łuku aż trzy znaki pokryją się ze sobą, co oznacza, że będzie to de facto lepsze połączenie, bo utworzy krótszą sekwencję.

Z powyższej obserwacji wynika, że należy zmodyfikować funkcję nadawania wartości łukom, aby uwzględniać relatywną 
liczbę powtarzających się znaków. W ten sposób tworzymy algorytm, który radzi sobie lepiej z instancjami zawierającymi błędy negatywne. Pogrubione zostały fragmenty, które uległy modyfikacji.

Algorytm 4
\begin{enumerate}
 \item Wczytaj zbiór wszystkich słów tworzących spektrum do słownika \texttt{OLIGS}
 \item Utwórz graf skierowany $G$. Niech zbiór wierzchołków \texttt{V} stanowią słowa w słowniku \texttt{OLIGS}. 
       Dla każdej pary wierzchołków $(i,j)$ dodaj łuk od wierzchołka $i$ do $j$ o koszcie równym 
       {\bf różnicy długości słowa $i$ do odległości słowa w wierzchołku $j$ względem słowa w wierzchołku $i$}
 \item Znajdź {\bf najwyższy koszt łuku \texttt{MAX}} w grafie $G$
 \item Utwórz graf skierowany $G'$, który jest podgrafem grafu $G$ i zawiera tylko łuki o koszcie $MAX$
 \item Usuń cykle w grafie $G'$ tworząc graf $G'_{DAG}$ korzystając z Algorytmu 1
 \item Uszereguj topologicznie $G'_{DAG}$
 \item Wybierz najdłuższą możliwą ścieżkę \texttt{PATH} w grafie $G'_{DAG}$
 \item Usuń ze słownika OLIGS słowa wchodzące w skład ścieżki \texttt{PATH}, w zamian dodając słowo utworzone przez połączenie słów zgodnie z kolejnością w ścieżce \texttt{PATH}
 \item Jeśli słownik OLIGS zawiera dokładnie jedno słowo, to zakończ, w przeciwnym wypadku idź do punktu 2
\end{enumerate}

Powyższy algorytm posiada identyczną złożoność jak algorytm ogólny, działa jednak znacznie lepiej.

\end{document}
